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In the second part of this paper in micro canonical ensemble the new numerical approach for 
consideration of quantum dynamics and calculations of the average values of quantum operators and 
time correlation functions in the Wigner representation of quantum statistical mechanics has been 
developed. The time correlation functions have been presented in the form of the integral of the 
Weyl's symbol of considered operators and the Fourier transform of the product of matrix elements 
of the dynamic propagators. For the last function the integral Wigner- Liouville's type equation has 
been derived. The initial condition for this equation has been obtained in the form of the Fourier 
transform of the Wiener path integral representation of the matrix elements of the propagators at 
initial time. The numerical procedure for solving this equation combining both molecular dynamics 
and Monte Carlo methods has been developed. 

The numerical results have been obtained for series of the average values of quantum operators 
as well as for the time correlation function characterizing the energy level structure, the momentum 
flow of tunneling particles at barrier crossing and the absorption spectra of electron in potential 
well. The developed quantum dynamics method was tested by comparison of numerical results 
with analytical estimations. Tunneling transitions and the effect of the quasi stationary state has 
been considered as the reason of the peculiarities in behaviour of the time correlation functions and 
position and momentum dispersions. 

Possibility of applying the developed approach to the theory of classical wave propagation in 
random media have been also considered. For classical waves some results have been obtained for 
Gaussian beam propagation in 2D and 3D waveguides. 
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I. INTRODUCTION 



In canonical ensemble considered in the first part of this paper |jj the numerical studies of the exponentially small 
tunneling effects is very difficult due to the temperature averaging. To overcome this difficulty the more delicate 
approach in the micro- canonical ensemble with a fixed initial energy has been introduced by the inverse Laplace 
transformation of the spectral density on the inverse temperature variable. The related integral Wigner- Liouville's 
type equation have been obtained and the numerical approach combining both molecular dynamics and Monte Carlo 
methods for solving this equation has been developed. 

The time correlation function characterizing the energy level structure, the momentum flow of tunneling particles 
at barrier crossing and the absorption spectra of electron in potential well have been calculated. Tunneling transitions 
and the effect of the quasi stationary state has been considered as the reason of the peculiarities in behaviour of the 
time correlation functions and position and momentum dispersions. 

Possibility of applying the developed approach to the theory of classical wave propagation in random media have 
been also considered. For classical waves some results have been obtained for Gaussian beam propagation in 2D and 
3D waveguides. 



II. ELECTRON IN DOUBLE WELL POTENTIAL 



In the second part of this paper we've considered dynamics of quantum electron in a deep symmetric (with respect 
zero of x-axis q x — 0) double well potential: 



VoU (\q\) = V {Vo'/Vb * exp (\q\ 2 /a 2 ) - exp (\q\ 2 /a 2 ) } 



where F ' < V and a > a (VJ/V = 0.12). 

To analyze the tunneling effects, absorption spectra and electron energy levels we have considered the Fourier 
transform of the time correlation functions characterizing the quantum particle momentum flow through the barrier 



k (a;, E) = J °° exp {—iwt — et) Cf-q (t, E) dt = 
Z~ x J °° exp {-iut - et) Tr (Pexp (iHt/h) fyexp l-iHt/h) 5 (e - H]) dt = 



(1) 



where e — > 0, 6 \ E — H^j is the initial density matrix , Z = Tr (s (e — Hj) > E i* are eigenvalues and eigenfunctions 

of the Hamiltonian of the system. Due to delta- functions the function k(u,E) should have peaks in discrete part 
of spectrum on the frequencies lo equal to the different combinations (E^ ~ E v ) /h. Here Cft/ (t, E) is defined in the 
Wigner representation by expression: 



Cf v (t, E) — 

J J d Pl d qi dp 2 dq 2 \ (F (pf , gf ) r, ((g) + F (pf , qi ) V (qf )) x 
W (p 1 ,q 1 ;p 2 ,q2;t;E) 



(2) 



where Weyl's symbols of operators F and r) are: F (p, q) = [p5 (q) + 5 (q)p\, r\ is step function that projects onto 
q x > half- space. 

Introduced in the first part of this paper jl) the spectral densities in canonical W (pi, qi',P2, 92; t; i%0) and micro- 
canonical W (pi, 5i;j>2) Q2) t; E) ensembles are connected according to the d|) by the Laplace transformation: 



W (p 1 ,q 1 ;p2,q2]t;ih(3) = J dEexp (-0E) W (pi, qi;p2, Q2; t; E) = 

J dE exp {—(3E) — J duo exp {iut — et) exp 11 



k(u,E) 



where 



W (pi,gi;p2,»; t; E) - Z' 1 J J d^dfr exp (i^f) exp (i^-) 



x (9i + % 



exp 



<12 



<12 



exp (-iHt/T^j 5(e - i?) 
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We suppose that function W (pi, <?i;P2, 12, t; E) provides more information about quantum effects and quantum 
dynamics of the system. 

The functions W (pi, qi',P2, 12\ t] E) as can be easily proved are the solutions of the same linear integral equation 
as have been obtained in the first part of this paper |lj . 

W(pi,qi;p 2 ,q 2 ;t;E) = W{p ,qo;Po,qo\F) + J dr J dsdr)W(p T - s,q T ;p T - r],q T ;r;E)^(s,q T ;r],q T ) (3) 

where j(s,q T ;r),q T ) = \{u (s,q T ) 6(rj) - uj (mq T ) S(s)}, uj{s,q) is 



w(s,q) 



(2 rt ).»/* v <«-«'» ,ta ( 2 r) + 'W * 



S(s) is the Dirac delta function, {q T (r;pi, qi, £), p T (T\px, <7i, t)} and {q T (r;p2, 92, t), p T (r; P2, <?2, t)} are pair of classical 
dynamical pq- trajectories for 'positive' and 'negative' time direction and initial condition taken at r = t: 

dp/dr = ^F(q T (r));q t (t;p 1 ,q 1 ,t) = q 1 

dq/dr = ^—p T (t) ; p t (t; p x , qi, t) = p x 
2m 

dp/dr = --F(q T (r));q t (t;p 2 ,q 2 ,t) = q 2 (4) 

dq/dr = -—p T (T);p t (t;p 2 ,q 2 ,t) = p 2 
2m 

The initial conditions W (pi, qv,p 2 , q 2 ] E) are connected with initial functions W (pi, (?i;p2, q2\ ihfS) in analogous 
way: 



W (pi,qi;p2,q 2 ;ih/3) = J dEexp (-(3E) W (pi, qv,p 2 , qr, E) 

f 1 f°° 

/ dEexp (—/3E) — / dujexp 

J 27T 



k {uj, E) 

For checking the basic ideas of this approach we have used the classical approximation of the initial condition of 
spectral density W (pi, qi;p 2 ,q2', E) in the following form: 

W {p 1 ,q 1 ;p2,q2-,E) « 5 (E - H (px,qi))rj{-qf) S (pi -p 2 )<5(9i - q 2 ) (5) 

where H (pi ,q x ) is the classical Hamiltonian of the system. The initial spectral density is non zero only in the half - 
space of negative part of the x - axis. 

Let us rewrite the integral equation and the iteration form of its solution in symbolic form: W t = W t + K l T W T 
and 

W* = W l + K* Ti W Tl + K^K^W^ + K^K^K^W^ + ... (6) 

Here W 1 and W Tl is the quantum initial density evolving classically in intervals [0,t] and [0, n], while K^l +1 are 
operators, which describe propagation between times 74 and r<+i. The time correlation functions are the linear 
functionals of the spectral density: 

C Fv (t, E) == J J dp 1 dq 1 dp 2 dq2\ (F (pf , qf) V (gf ) + F (pf , g£) V (gf )) x 

W(pi,5i;p2,<&;t;B) = (7) 
+ (</>|i^ Tl ) + (^|i^ 2 ^^) + fa|l£*,iq*jq>W™) + ... 

where <l>(pi,qi;p 2 ,q 2 ) = 5 (F (pf, gf ) 77 (gf) + F (pf, gf ) 77 (gf )), brackets (|) for functions and W"(po, 9o;po, <7o; or 
K*.K^ i l ...K^W T1 mean the integration over the phase spaces {pi,gi;p2,g2}- 

Introduced in the first part of this paper jjj the recurrent relations for the pieces of dynamic trajectories according 

to a 



3 



=3+1 =7+1 \ 

3 - + 1 o. ^J+i 



2 ,/ 
p 



-/ 4+1 4 + 1 A W 

^j+i „. sj+1 



allow to obtain the explicit expression of the terms of series (^J) and to analyze its mathematical structure. So for 
example the third term (J = 3) can be written as: 



dr 2 \ 
i) Jo 



X7(«2, 925 V2, q2)l{si,ql; Vi,q~iW(pl, ql;pl, ql\E) 
Note that average values of quantum operators A (t) can be formally presented in the form analogous to (Q) : 

A (t, E) = Z~ l Tr (exp (iHt/fij iexp (-iHt/h\ 5 (e — Hj 



xVF (p 1 ,q 1 ;p 2 ,q 2 ;t- 1 E) 



III. WIGNER APPROACH IN THE THEORY OF CLASSICAL WAVE PROPAGATION IN RANDOM 

MEDIA 

In the case of the scale of inhomogeneities is large in comparison with the wave length it is possible to neglect 
the backscattering and depolarization of waves and to describe the classical wave propagation by the parabolic wave 
equation Q : 

du 1 

i— = --An — e (z, R) u; u (z , R 1 ui) = u° (R, ui) (9) 

where fc _1 is taken as a unit of length, k — ui (e) 3 /c, u> is the wave frequency, z is the initial direction of wave 

propagation, A = {-§^) 2 + C^jj) > ^ ( 2 ; R) = ( £ ~ £ ) l^i £ ' IS ^ ne dielectric permittivity, the overbar denotes averaging 

over the medium fluctuations, u(z,R,oj) is the slowly varying complex amplitude of the wave field [Q , u° is the 
'initial condition' for u at z = zo and R = (x, y) is a radius- vector in transverse plane to the initial direction of the 
wave propagation (axis z ). The Wigner- Liouville distribution function is defined by 

f(z,P,Q,u) = J 2^ r Jp(z,Q-l,Q + l,u) exp(zPO^; 
p (z, Q — |, Q + |, u?j = u (z, R\, lj) u* (z, i?2, w) = u (z, Q — §, to) u* (z, Q + |, (J) 

Taking the time derivatives of the Wigner- Liouville distribution function and using equation (^) it is possible to 
obtain the integral equation for / (z, P, Q, ui) in the form: 



/ 0, P, Q, w ) = f° (zq, Po,Qo, w) + Jo dr J dSf (r, P T - S, Q T , ui) w (r, S, Q T ) , 
w (r, S, Q) = -j^tyrK I d Q' 1 ( T > Q - Q') sin (^^f) + F r (Q) ^ 

where t G [0, z\ . The dynamic trajectories \Q t [t\ P, Q, z), P t (t; P, Q, z)} with 'initial conditions' {P, Q} at r = z 
are defined by the following equations: 

dP/dT = F T {Q T {T))- Q z (z;P,Q,z) = Q 

dQ/dT = P T {T)/m; P z (z;P,Q,z)=P (10) 
F T (Q)=d(e(r,Q))/dQ 
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So in parabolic approximation the variable z can be formally considered as the "time variable" in the Schrodinger 
like equation and the developed approach can be used for calculations of the average intensity, wave fields moments, 
the scintillation index and correlation functions characterizing the classical wave scattering and propagation in random 
media. We'll restrict our interests with linear functionals of Wigner function, representing the mentioned above values 



(11) 



where brackets (|) mean integration of functions <f> an d f (z,P,Q,w) or K^Kr J j _ 1 ---K^ f Tl over phase space (P,Q). 
This represents also the averaged value of operators (f> I P, Q ) in the Wigner representation 



. (z, w) = — / dPdQ^ (P, Q) f (z, P, Q, u) 



,(P.()) = j dZexp(iPd/h)(Q-^ 



P,Q 



Q 



where ( Q — | 



(p> Qj Q + I y is a matrix element of operator 4> [§ ■ 



(12) 



IV. QUANTUM AND WAVE DYNAMICS 



The possibility to convert series like (J7]) and (^Tj) into the form convenient for probabilistic interpretation allow us 
to develop the Monte Carlo method for its calculation || , Q . Let us note only that the ergodic hypothesis allow us 
to perform the integration on variables {p%, qi',P2, 12} for averaging in micro canonical ensemble according to (j5|) , (Q) 
by using the classic molecular dynamics method. As has been shown in Q the transition from the integration over 
the phase- space at the fixed total energy (due to (||)) to the time averaging along the classical dynamic trajectory 
requires the correction factor \\/H (pi, qi)\, which in our case should be added to the weight function £1 of the p - 
trajectories j^], ||. Here multidimensional vector \/H (pi,qi) has the following components 

{dH/dp\, dH/dp\,dH/dq\, OH/dqX, 
dH/dp^, ...,dH/dp v N ,dH/dq l N , ...,dH/dq v N } 

for the v - dimensional space, while the dynamic evolution has been realized according to equations (S). 



V. NUMERICAL RESULTS 



A. Electron spectra 



The Fig. |l| shows the time correlation function CFri(t,E) for E/Vq — —.92. From physical point of view this 
function characterizes the momentum flow of quantum particles tunneling through the potential barrier. First of all 
it is necessary to note that at small (less than 20 th/Vo) and large times (more than 20 th/Vo) the time correlation 
function has two different frequencies of oscillation. Physically this behaviour can be explained by the shape of our 
double well, which have two space scales and two energy scales. The left and right shallow narrow wells separated by 
barrier have the space scale of the order a and the depth equal to 0.12 * Vb- The space scale of the large well containing 
both shallow narrow wells and the barrier is a, while the depth is Vq. These scales satisfies as we mentioned before 
to the following inequalities: 0.12 = Vq/Vo <C 1 and a > a. 

The main contribution to the time correlation function at initial time (less 20) comes from the quantum trajectories 
with virtual energy close to the height of the barrier, when they path trough the top of the barrier. Due to the 
momentum jumps these trajectories can penetrate in both shallow wells and can be also trapped there. So the high 
frequency oscillations can be connected with multiple reflections and tunneling transitions of these trajectories . It 
is necessary to stress that derivative of this correlation function with respect to time is the momentum- momentum 
time correlation function taken at the top of the barrier. So from physical point of view these oscillation results from 
changing of the direction of the main momentum flow of quantum particle at multiple tunneling transitions. 
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The same functions obtained in approximation of the classical trajectories (the first term of the iteration series) 
are identically equal to zero. Indeed at the averaging in the micro- canonical ensemble the initial p, q data for our 
classical and quantum trajectories have been taken from the classical trajectory moving in the left part of our double 
well. That means that the initial energy of our trajectories was lower than the height of the barrier and the classical 
trajectories without momentum jumps were unable to leave the left part of our double well and consequently to give 
contribution to the considered time correlation function. 

At the large time (more than 20 th/Vo) the main contribution to the time correlation function considered comes 
from the trajectories with the virtual energy far (much more) from the height of the barrier. So these trajectories can 
move anywhere in the large well and that is the reason of low frequency oscillations of the time correlation function. 

The next Fig. || and Fig. || present the squared amplitude of the Fourier transform of the time correlation function 
\k {o-i,E) | 2 versus the frequency Hlj/Vq. Analytical semi classical estimations of the positions of the sharp peaks of 
the function k (oj, E) are presented by stars 1. Calculated positions of the sharp peaks of the \k (uj, E) | 2 on the Fig. 
|2| and Fig. |3| are in a agreement with analytical estimations. 

B. Position and momentum dispersions at quantum tunnelling 

The next Fig. [| and Fig. || present results for position dispersions obtained for classical and quantum trajectories 
respectively versus dimensionless time tV^/h for ID case. As we mentioned before the initial energy of the our 
trajectories was lower than the height of the barrier and calculations of position dispersion allowing for only the first 
term of the iteration series were implemented by making use of only the classical trajectories, which were unable to 
leave the left shallow well. These results for position dispersion are presented on the Fig. The squared root of the 
position dispersion gives the estimation of the oscillation amplitude (of order of 0.05/c -1 ) of these trajectories and 
allow to estimate the characteristic size of the available space in the left potential well for these trajectories. Note 
that the classical dynamics gives non damping oscillations of position dispersion. 

The quantum trajectories give the qualitatively different behaviour of position dispersion. These results are pre- 
sented on Fig. |^. The virtual energy of quantum trajectories can be larger than the height of the barrier, so these 
trajectories can move in the right part of the double well, where the virtual energy of these trajectories can become 
smaller than the height of the barrier and these trajectories can be trapped there. At any case the quantum trajec- 
tories can travel anywhere in the double well and even leave it. So the position dispersion of quantum trajectories 
is much larger than the same value of classical trajectories. However position dispersion of the quantum trajectories 
has a very interesting peculiarity. At time tVo/h = 30 the position dispersion becomes practically equal to zero. It is 
interesting that this happens at time when time correlation function changes the frequency of its oscillations. Note 
that at this time the momentum dispersion presented on next Fig. ^ (curve 2) changes also its characteristic behaviour 
after sharp oscillation. One can see that the minima of position and momentum dispersions happens approximately 
at the same time. Then the momentum dispersion practically stabilizes. 

Our physical explanation of this peculiarity is connected not only with tunneling transitions of quantum particles 
but also with existing of a quasi stationary state at the top of the narrow barrier. This state may be the reason that 
both position and momentum dispersions at the same time (of about 30 tVo/ft, ) have minima. The time of life of this 
quasi stationary state may be estimated from the Fig. || as of order about htVo/H 

The next Fig. [| presents the momentum dispersion for classical and quantum trajectories. The momentum 
dispersion for classical trajectories oscillates with very small amplitude (curve 1). 

VI. CLASSICAL WAVE PROPAGATION IN 2D AND 3D WAVEGUIDES 

To test the developed stochastic dynamics approach the wave propagation along z-axis (z > 0) in 2D and 3D 
waveguides with realistic profile of refractive index has been investigated in parabolic approximation. Gaussian beam 
distribution was used as a initial condition at z = 0. Figures |?] and || present the 3D case numerical data for average 
position (x (z) , y (z)) and dispersion (3 = R (z) 2 — R (z) 2 of Gaussian beam vs distance z. Foci points along waveguide 
are indicated by the minimum values of (3. 

VII. CONCLUSION 

In the Wigner formulation of quantum statistical mechanics for canonical and micro canonical ensembles we have 
presented a new computational technique allowing quantum dynamics simulations for systems including subsystems of 



G 



quantum interacting particles and subsystems of classical heavy scatterers as well as the system of quantum particles 
in external potential field. The developed approach for quantum dynamics includes a sophisticated combination of 
well known molecular dynamics method and Monte Carlo technique. Numerical results have been presented for the 
time correlation function characterizing the energy level structure, the momentum flow of tunneling particles at barrier 
crossing and the absorption spectra of electron in potential well have been calculated. Tunneling transitions and the 
effect of the quasi stationary state has been considered as the reason of the peculiarities in behaviour of the time 
correlation functions and position and momentum dispersions. 

Possibility of applying the developed approach to the theory of classical wave propagation in random media have 
been also considered. For classical waves some results have been obtained for Gaussian beam propagation in 2D and 
3D waveguides. 
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FIG. 1. The time correlation function Cft, (t,E) versus tVo/Ti for quantum trajectories. (E/Vo = —0.92) 
FIG. 2. \k(u>,E) | 2 vz frequency hu>/Vo- 1- analitical estimations of the sharp peak positions; 2 - numerical results; 3- the 
error bar; a = 0.26, a = 2.23. 

FIG. 3. \k (u>, E) 2 vz frequency hu>/Vo- 1- analitical estimations of the sharp peak positions; 2 - numerical results; 3- the 
error bar; a = 0.26, a = 3.43. 

FIG. 4. Position dispersion: 1- classical trajectories; (E/Vq = —0.92) 

FIG. 5. Position dispersion: 1- quantum trajectories; (E/Vo = —0.92) 
FIG. 6. Momentum dispersion: 1- classical trajectories; 2- quantum trajectories; (E/Vo = —0.92) 

FIG. 7. Average position x(z) and y(z) vs z for 3D waveguide. Starting point at z = corresponds approximately to point 
(21,21), 'center' of waveguide is nearly the point (17,17). 

FIG. 8. 13 = R 2 ~(z) - (R{z)) 2 vs z (in discrete units ) for 3D waveguide (R =(x,y)). 
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